---
title: "Power Analysis"
author: "Kaylyn Jackson Schiff, Daniel Schiff, and Natalia Bueno"
date: "2020"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

#####Setup Chunk
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rmarkdown)
library(tidyverse)
library(estimatr)
options(scipen=999)
```


#####Set Parameters for Power Analysis
```{r establish parameters for power analysis}
alpha <- 0.05 # Standard significance level 
possible.n <- seq(from=500, to=3000, by=100) # The sample sizes we'll be considering 
sims <- 1000 # Number of simulated treatment allocations to conduct for each N 

powers <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
powers.diff <- rep(NA, length(possible.n)) # Empty object to collect lm simulation estimates 
```


####Simulate Results for Politican Response on Support Outcome (Liar's Dividend)
```{r loop of lms to simulate results}
#Loop to vary the number of subjects
for (j in 1:length(possible.n)){ N <- possible.n[j] # Pick the sample size

Y0 <- rnorm(n=N, mean=-.0722, sd=.924) # control potential outcome 
tau1 <- .1776 # Hypothesize treatment effect for IU
Y1 <- Y0 + tau1 # treatment potential outcome for IU
Y2 <- Y0 + 0 # treatment effect for fact-checking


significant.lm <- rep(NA, sims) # Empty object to count significant experiments 

#Inner loop to conduct experiments "sims" times over for each N
for (i in 1:sims){
  
  Z.sim <-  sample(0:2, N, replace = T, prob = c(1/3, 1/3, 1/3)) # Do a random assignment 
  Y.sim <- case_when(Z.sim == 0 ~ Y0, Z.sim == 1 ~ Y1, Z.sim == 2 ~ Y2) # Reveal outcomes according to assignment 

  fit.lm <- lm_robust(Y.sim ~ factor(Z.sim)) # Do analysis (Simple regression)
  
  p.lm <- summary(fit.lm)$coefficients[2,4] # Extract p-values  
  significant.lm[i] <- (p.lm <= alpha) # Determine significance according to p <= 0.05
}
  powers[j] <- mean(significant.lm) # store average success rate (power) for each N 
}

powers1 <- powers
```

####Plot Results
```{r plot results}
png("Figures/power_support_followup.png")
plot(possible.n, powers1, type = "b", main = "Power to Detect IU Effect on Support by Sample Size", 
     xlab = "Sample Size", ylab = "Power")
points(possible.n, powers1, 
       pch = 19, col = "gray", bg = "gray", cex = 1.3)
abline(h = .8, lty = 2)
#abline(v = 1500, lty = 4)
dev.off()
```

